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DECONVOLUTION TECHNIQUE EMPLOYING HERMITE FUNCTIONS 

STATEMENT OF GOVERNMENT INTEREST 
Under paragraph 1(a) of Executive Order 10096, the conditions under which this 
5 invention was made entitle the Government of the United States, as represented by the 
Secretary of the Army, to the entire right, title and interest therein of any patent granted 
thereon by the United States. This and related patents are available for licensing to 
qualified licensees. Please contact John Griffin at 703 428-6265 or Phillip Stewart at 601 
634-4113. 

10 

FIELD OF THE INVENTION 
The present invention relates generally to mathematical processes for removing 
distortion from raw data that presents in a Gaussian distribution. In particular, it describes 
an efficient improved de-convolution procedure for data processing adapted according to 
15 specific needs of the user. It applies to all such measurements including optical and 
electro-optical measurements. 

BACKGROUND 

Generally speaking, a measurement or imaging system introduces "distortion" 
20 into the "measurement" or image that is undesirable, i.e., it distorts the desired 
information. This occurs in all measuring systems. This undesirable distortion may be 
from a number of sources including the characteristics of the measuring system itself, the 
environment in which the measurable or image is embedded, the environment between 
the detector(s) of the measurement system and the item to be measured or imaged, the 
25 experience of the operator of the measuring system, etc. The resulting agglomeration of 
undesirable distortions that have been "added" to the actual data results from a 
"convolution" of one or more of these sources of undesirable information with the desired 
(actual or undistorted) data by the measuring or imaging system. 

The process for getting actual (undistorted) data from measured (distorted) data 
30 involves implementing a "deconvolution" that attempts to reduce or eliminate the 
undesirable distortion. To do this, deconvolution techniques and deconvolution filters 
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have been developed, some examples of which are represented in the following U. S. 
patents and published applications. 

A number of U.S. patents involve some form of deconvolution technique or filter 
for which the methods of the present invention may provide an improvement. These 
5 include the following examples. 

U. S. patent no. 5,185,805, Tuned Deconvolution Digital Filter for Elimination of 
Loudspeaker Output Blurring, to Chiang, February 9, 1993, is incorporated herein by 
reference. The response of a speaker to a known analog signal is sampled at or above the 
Nyquist rate and the result is used to construct a linear, time invariant deconvolution filter 
10 that compacts the blurred output signal back into its original (unblurred) form. The 
frequency response is then improved by a fine-tuning process using Lagrange's Method 
of Multipliers. 

U. S. patent no. 5,400,265, Procedure for Enhancing Resolution of Spectral 
Information, to Kauppinen, March 21, 1995, is incorporated herein by reference. The 
15 Fourier Self-Deconvolution (FSD) method is used to enhance resolution of spectral data. 
The Maximum Entropy Method (MEM) is then used to predict filter error coefficients. 
Using data from the FSD and the MEM and the Linear Prediction (LP) method, data are 
predicted beyond the area provided by the FSD to maximize spectral line narrowing with 
minimal distortion. 

20 U. S. patent no. 5,400,299, Seismic Vibrator Signature Deconvolution, to 

Trantham, March 21, 1995, is incorporated herein by reference. A deterministic signature 
deconvolution of the data compresses the impulse response of the data resulting in 
sharper seismic images than possible using cross correlation techniques. 

U. S. patent no. 5,748,491, Deconvolution Method for the Analysis of Data 
25 Resulting from Analytical Separation Processes, to Allison et al., May 5, 1998, is 
incorporated herein by reference. A signal representing multiple partially separated 
sample zones is measured, a Point Spread Function (PSF) determined, and a Fourier 
transform of the data performed. A noise component of the signal is determined and a 
value for a resulting signal is determined using a specific filter. The inverse Fourier 
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Transform is then taken to get the desired result. Although a Gaussian PSF is preferred 
for use in the process, a method is provided for surveying a number of PSF possibilities 
serially. 

U. S. patent no. 5,761,345, Method of Discrete Orthogonal Basis Restoration, to 
5 Moody, June 2, 1998, is incorporated herein by reference. The method applies to linear 
data and involves estimating a signal-to-noise ratio (SNR) for a restored signal or image. 
A set of orthogonal basis functions is then selected to provide a stable inverse solution 
based on the estimated SNR. Time or spatially varying distortions in the restored system 
are removed and an appropriate inverse solution vector obtained. 

10 U. S. patent no. 5,862,269, Apparatus and Method for Rapidly Convergent 

Parallel Processed Deconvolution, to Cohen et al., January 19, 1999, is incorporated 
herein by reference. Provided is a stand-alone image processing system and method for 
deconvolution or an initiator for rapid convergence of conventional deconvolution 
methods. The method improves processing speed by relaxing the sequential requirement 

15 of the CLEAN method. The number of iterations is reduced by fractional removal of 
noise for multiple image features during the processing of a single subtractive iteration. 

U. S. patent application publication no. 2002/0156821 Al, Signal Processing 
Using the Self Deconvolving Data Reconstruction Algorithm, by Caron, October 24, 
2002, is incorporated herein by reference. A filter function is extracted from degraded 
20 data by mathematical operations applying a power law and smoothing function in 
frequency space. The filter function is used to restore the degraded content through any 
deconvolution algorithm without prior knowledge of the detection system, i.e., blind 
deconvolution. 

U. S. patent application publication no. 2002/0136133 Al, Superresolution in 
25 Periodic Data Storage Media, by Kraemer et al., September 26, 2002, is incorporated 
herein by reference. A method termed Matrix-Method Deconvolution (MMD) 
compensates for the effects of an optical addressing system's PSF. Prior knowledge of the 
PSF and inter-memory center spacing is used to acquire binary data stored physically in a 
periodic storage medium. 
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Optical and electro-optical devices cause intrinsic distortion of the image or 
spectral characteristic they sense, measure or record. This distortion is separate from any 
measurement noise or artifact noise that may be present. The mathematical relationship 
between the undistorted and "intrinsically distorted" images or spectra is described by an 
5 equation incorporating a convolution integral. Blass, W. E. and G. W. Halsey, 
Deconvolution of Absorption Spectra, Academic Press, New York, 1981. Jansson, P. A., 
Ed., Deconvolution of Images and Spectra, Academic Press, New York, 1997. 

The measured (distorted) image equals the convolution integral of the product of 
the actual (undistorted) image and the point spread function (PSF) that stochastically 

10 characterizes the device. (Consider a very small point of light. If the visual system had 
perfect optics the image of this point on a lens would be identical to the original point of 
light. Referring to Fig. 1, if the relative intensity of this point of light were plotted as a 
function of distance on the lens, such a plot would look like the dashed vertical line. 
However, the optics are not perfect so the relative intensity of the point of light may be 

15 distributed across the lens as shown by the curved line. This curve is the PSF.) 

Conventionally deriving an undistorted image or spectra from an intrinsically 
distorted one may be attempted initially by a stand-alone process termed classical 
deconvolution. Unlike the convolution integral, deconvolution is not a well-defined 
mathematical operation, i.e., only the desired goal is well defined. 

20 Because of the disadvantages of the deconvolution process, some have offered 

solutions to completely avoid it. See, for example, U. S. patent no. 3,934,485, Method 
and Apparatus for Pulse Echo Imaging, to Beretsky et al., January 27, 1976. 

Some authors use the term "deconvolution" for what is actually a combination of 
classical (intrinsic) deconvolution and one of several special iteration techniques. Blass 
25 (1981). Jansson (1997). Mou-yan, Z. et al., On the Computational Model of a Kind of 
Deconvolution Problem, IEEE Trans. Image Process , 4, No. 10, pp. 1464-1476, Oct 
1995. This is a long-standing problem in deconvolution research, i.e., extracting an 
undistorted image from solutions to incompatible systems of linear equations derived 
from "discretizing" (not merely digitizing) convolution integrals. This discretizing yields 
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a system of linear equations that relate the actual (undistorted) and captured (distorted) 
images in deconvolution. Twomey, S., Introduction to the Mathematics of Inversion in 
Remote Sensing and Indirect Measurements, Elsevier Sci. Pub., pp. 33-35, 1977. 

Many classical deconvolution techniques, when used alone, produce large 
5 unacceptable errors. To reduce these errors, the aforementioned linear iteration 
techniques are used with intrinsic deconvolution. Used together, the combination, after 
many iterations, may produce acceptable error levels. Use of these combinations also 
may increase signal-to-noise ratios (SNRs). Blass (1981). Jansson (1997). Mou-yan et al. 
(1995). 

A recent patent details use of the Hermite Functions to improve performance of a 
sonar system. U. S. patent no. 6,466,515 Bl, Power-Efficient Sonar System Employing a 
Waveform and Processing Method for Improved Range Resolution at High Doppler 
Sensitivity, to Alsup et al., October 15, 2002, is incorporated herein by reference. A sonar 
system incorporates a comb-like waveform by modulating its lines according to a set of 
Hermite functions defining a Hermite Function Space (HFS). Doppler sensitivity akin to 
CW systems is realized by applying to the HFS signals, the deconvolution method of the 
'515 patent. Although Hermite Functions are used, they are not used to simplify the 
deconvolution of received data. 

Thus, although a number of "combined" deconvolution techniques are known, 
20 difficulties arise in applying many of them. The present invention provides a solution that 
overcomes such difficulties. 

SUMMARY 

A unique deconvolution procedure adapts to generate unique deconvolution 
algorithms for each different requirement. This is accomplished by first solving the 
25 general convolution equation exactly, in closed analytical form. Mou-Yan, Z. et al., New 
Algorithms of Two-Dimensional Blind Deconvolution, Opt. Eng., No. 10, pp. 2945-2956, 
Oct. 1995. Next, the results are transformed, leading to a fundamentally new linear 
relationship between the actual (undistorted) and the captured (distorted) images. 
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To do this, a preferred embodiment of the present invention uses the Hermite 
functions and the Fourier-Hermite series to represent the two classes (distorted and 
undistorted) of images. Using this procedure leads to the exact solution of the 
convolution integral. 

5 A preferred embodiment of the present invention provides a way around the need 

for solving incompatible systems of linear equations derived from "discretizing" 
convolution integrals. No numerical evaluation of the convolution integral is required in 
the present invention. 

A preferred embodiment of the present invention is executed by exploiting a 

10 remarkable mathematical coincidence. The mathematical coincidence is that the most 
common PSF used to characterize a device is a Gaussian function (see Fig. 1). The 
Gaussian function is also a Fourier-Hermite function of zero order. 

Exploitation of this fact permits exact analytic evaluation of the convolution 
integral for most optical imaging systems and some non-optical systems. The 

15 mathematics is processed in exact closed form, followed by analytic deconvolution. 
Further, the present invention addresses an additional concern. When employing 
conventional "combined deconvolution," image distortion may also result because of the 
numerical division of very noisy data represented by PSF "points" having values less than 
one. Division in the present invention is by a new function of the PSF parameters 

20 yielding values generally greater than one. 

If desired, the procedure may be used with any of the iterative techniques to 
further refine the image. Finally, if used with an iterative technique, fewer iterations will 
be needed because of the inherently more accurate starting point provided by the 
procedure. 

25 Specifically provided is a process for deconvolving data collected by a device the 

point spread function response of which follows a Gaussian distribution. The process 
accurately estimates actual parameters derivable from this "captured" data. It comprises: 
forming a mathematical relationship having a first mathematical equivalent of the 
captured data on one side of an equality and a second mathematical equivalent of 
30 actual parameters on the other side of the equality; 
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selecting an order, w, of a Hermite function for modifying this mathematical 
relationship; 

modifying the mathematical relationship to form Hermite functions therein to 
permit identification of like items, each having a coefficient, on each side of the 
5 equality, the coefficients associated with the actual parameters being unknowns; 

developing a set of linear equations from this mathematical relationship that relate 
to the coefficients of the like items; and 

solving this set of linear equations for the unknown coefficients to produce an 

exact solution to the convolution; and 
10 deconvolving this set of linear equations to determine the unknown coefficients. 

The process may also initiate any conventional iterative deconvolution techniques 
to further refine the data, such that fewer iterations will be needed because of provision of 
an accurate starting point. 

Specifically, a preferred embodiment of the present invention is a process for 
15 deconvolving output from detectors whose point spread function follows a Gaussian 
distribution. The process accurately estimates environments derivable from the detectors' 
output. The process comprises: 

forming a mathematical relationship having the mathematical equivalent of the 

output on one side of an equality and the mathematical equivalent of the 
20 environments on the other side; 

selecting an order, w, of a Hermite function for modifying this equality, such that 

m also determines the number of terms for the representation of the environments; 

modifying the mathematical relationship to form Hermite functions therein, such 

that forming the Hermite functions permits identification of like items, each 
25 having a coefficient, on each side of the equality; 

developing a set of linear equations from the equality that relate to the coefficient 

of each of the like items; and 

solving the set of linear equations for the unknown coefficients of the like items to 
produce an exact solution to the convolution, in turn, permitting deconvolution 
30 using simple linear relationships. 
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Also, this embodiment may be used to initiate any conventional iterative 
nonlinear deconvolution techniques to further refine the output such that fewer iterations 
are needed because of the provided accurate starting point. 

In yet another preferred embodiment of the present invention is provided a 
process yielding accurate representations of actual image data by deconvolving image 
data collected by optical detectors whose point spread function response follows a 
Gaussian distribution. The process comprises: 

establishing a mathematical relationship as a general optical convolution integral 

equating the actual image data to the collected image data; 

selecting a Fourier-Hermite function; 

employing a generating function for Hermite polynomials for expanding a 
Fourier-Hermite series of the Fourier-Hermite function to establish a linear 
mathematical relationship between the actual and the collected image data; 
selecting an order, m, of the Fourier-Hermite polynomial such that m also 
determines the number of terms for the representation of the images; 
expanding the actual image data in a Fourier-Hermite form with unknown (to-be- 
determined) coefficients by employing a series of special transformations to 
convert the side of the mathematical relationship representing the actual image to 
a Fourier-Hermite series; 

expanding the collected image data in a Fourier-Hermite form; 

equating coefficients of like terms on each side of the mathematical relationship 

to relate the coefficients of the actual and collected images; 

selecting a linear convolution algorithm from which a general convolutional 

equation is derived; 

solving the resultant general convolution equation exactly, in closed analytical 
form, as a linear convolution equation, such that using this process to evaluate the 
general convolution integral yields a form proportional to a Gaussian function 
times a power series that is defined with a finite number of terms, m, and 
incorporates the unknown coefficients of the actual image data as presented in a 
Hermite function, such that a solution in closed analytic form provides a 
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satisfactory solution to the general convolution integral without approximation; 
and 

performing an analytic deconvolution of the resultant linear convolution 

algorithm by employing classical deconvolution algorithms solely, such that the 
5 analytic deconvolution yields a solution having acceptable error levels. 

As with other embodiments, this process may be used to initiate any conventional 
iterative nonlinear deconvolution techniques to further refine the accurate representations 
and fewer iterations are needed because of the accurate starting point provided. 

A preferred embodiment of the present invention for deconvolving captured 
10 (distorted) optical, or other, measurements has the following advantages: 

it completely circumvents the long-standing problem of having to use solutions 

from incompatible systems of linear equations derived from discretizing 

convolution integrals; 

it eliminates concern about discernible image distortion introduced by the imaging 
15 process itself in addition to that from the measuring device, e.g., a detector array; 

it automatically provides the unique algorithm(s) needed for the specific 
requirements of the user; 

it solves both the convolution and the deconvolution algorithms exactly through a 
system of related linear equations, not as an estimate based on an iterative 
20 process; 

it may be used with existing iterative processes, reducing the number of iterations 
required while further refining an image, if required; and 

it capitalizes on the unique characteristics of Hermite functions, enabling 

calculations in the frequency (wave number) domain as easily as in the time 
25 domain. 

Applications in the imaging field include: 

straightforward calculation of the actual (undistorted) measurement from the 

collection system's (distorted) measurement; 

removal of distortions from image data; and 
30 assessment of the accuracy of a given convolution. 
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Further advantages of the present invention will be apparent from the description 
below with reference to the accompanying drawings, in which like numbers indicate like 
elements. 

5 BRIEF DESCRIPTION OF THE DRAWINGS 

Figure 1 is a representation of a point spread function (PSF) compared to a theoretically 
perfect representation of a single point of light as may be captured by a lens. 

10 Figure 2 is a flow chart representing general steps used in implementing a preferred 
embodiment of the present invention. 

DETAILED DESCRIPTION 
There is concern about some of the numerical division operations that are a part 

15 of deconvolution algorithms. Concerns center on the division of a value associated with a 
signal datum that is mostly noise by the actual value of the PSF when the PSF is a very 
small percentage of its peak value, such as in the tails of a PSF represented by a Gaussian 
distribution. Of course, the peak value of the PSF, by definition, is one. It can be seen that 
such operations, i.e., division by a number much smaller than one, will magnify the 

20 contribution of the noisy datum. This produces significant distortion as a result of 
processing alone. This is then combined with the physical distortion introduced by the 
collection device itself. 

The Hermite- function method, proposed as a preferred embodiment of the present 
invention, prevents, or at least significantly reduces the chances of, these small values 

25 being used as a divisor. This algorithm provides for dividing by a function of the 
parameters of the PSF, as opposed to the actual values of the PSF. Thus, proposed is a 
function of the PSF \f(PSF) f independent of the point at which the PSF is being used. For 
most devices, this function is rarely much less than one, leading to little chance of 
magnifying a noisy datum. 
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The actual (undistorted) measurement is expanded in a Fourier-Hermite form with 
to-be-determined coefficients. This leads to the need to evaluate convolution integrals of 
the form: 



displaced PSF is represented by e 2 , 

the Hermite Function is represented by e 2 H m (x) , and 
H m (x) is the Hermite polynomial of order m. 

By way of introduction to the development of the theory behind the present 
invention, the "convolutional" relationship between the true and distorted measurements 
(or images) is given by: 

The recorded (distorted) image may be expressed as a convolution integral of the 

form: 



Suppose the convolutional equation relating the distorted and undistorted images 
could be solved exactly. Further, suppose the two sides of the equation can be 
manipulated to be of the same form. Thus, equating coefficients of like terms directly 
relates the parameters of the two images (distorted and actual), ultimately determining the 
desired algorithm. 

An expression for an actual (undistorted) image, I(x)> in terms of a Hermite 
function is: 




(i) 



where: 



PSF is represented by: Pe 1 , 0 < P < 1 , 




(2) 



(3) 



where H m (y) is the Hermite polynomial of order, m; and 



y = a 2 x 



(4) 
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where: 

x = the spatial variable and 

a = an arbitrary parameter 
Pauling, L. and E. Bright Wilson, Introduction to Quantum Mechanics, p. 80, McGraw- 
5 Hill, 1935. The Pauling reference includes a normalization factor that is a function of a. 
Since the factor is a constant, it has been absorbed in 7 m . 

The general convolutional relationship between the actual (undistorted) and 
captured (distorted) measurements (images) is given by: 

D(z)=jp(z-x)I(x)dx (5) 

10 where: 

D(z) = captured measurement; 

I(x) = actual value of the measurement; 

p(x) = point spread function (PSF) characterizing the measuring device, 
such as an optical or electro-optical imaging system 

15 

In most cases, p(x) may be described by: 

M) 2 

p(x) = Pe la (6) 

where: 

P is (2nb) t; \ the multiplicative "amplitude" of the Gaussian; 
20 a is a constant related to a measure of the "width" of the Gaussian; 

b is a constant that is an inverse measure of the "amplitude" of the Gaussian; 
and 

d is the displacement of Gaussian peak from the origin, normally zero 

25 Thus, the resulting form of P is in terms of a and b (d normally zero). When a = b, the 
integral of p(x) over all x is unity. P can never be greater than unity and for simplicity it is 
chosen as unity. 

When the PSF can be modeled by a Gaussian curve written as 

p(x) = Pe 2a (7) 

30 where: 
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a = a constant, 
b = a constant 
Expanding the unknown undistorted measurement (actual image), I(x), and the measured 
5 (distorted image), D(R), in Fourier-Hermite representations: 

ax 2 oo J_ 

/(*) = *' 2 £/A(«r 5 ) (8) 

and 

D{R) = f,D m e~H m (R) (9) 

10 7 m represents unknown coefficients that when solved for represent the solution of 

the deconvolutibn problem. D m represents known coefficients needed to represent D(R). 
Also, a is and arbitrary parameter that is later constrained by the requirement of physical 
"readability" of the description. 

D(R) and I(x) are related through the definitive equation containing the 

1 5 convolutional integral 

D(z) = £ p(z-x)I(x)dx (10) 

Integrating the right hand side (RHS) of Eqn. (10), a polynomial, I m (R) m , results. This 
polynomial must be transformed into the Hermite Polynomial form, J m H m (R), where J m 
is related to I m by a linear transformation the details of which depend on the order of the 
20 polynomial, m. 

The results produce identical structure on both sides of the equation, enabling the 
coefficients of like terms to be equated. A set of linear equations relating parameters of 
the distorted and undistorted images may now be written. Writing: 



(11) 



25 yields a system of linear equations: 



D m =a m (12) 
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10 



15 



20 



C = 



K 2bj 



that may be used to solve the convolution problem. 

Further, Eqn. (11) and (12) may be used to solve for J m and hence, I m , yielding the 
solution to the deconvolution problem. The common divisor, C, or a common multiplier, 
E, may be used, where E=l/C. Thus, 



E = 



26 
a ) 



so that 



J =ED 

m m 

Thus, the solution to the deconvolution problem can be expressed as 



(13) 



(14) 



(15) 



m=0 »i=0 

A process leading to a relationship equivalent to the integral of Eqn. (1) but able 
to be easily evaluated establishes a preferred embodiment of the present invention. The 
artifice of using a generating function for Hermite polynomials enables this. 

One generating function for Hermite polynomials is: 



T(x,t) = e- 2+2a =X 



ml 



(16) 



where / is a dummy variable. 

A usable integral may be created using the generating function as follows: 

('-*)' ± 

n = ]T(x,t)e 2 e 2 dx (17) 
Substituting the exponential definition in Eqn. (1) into Eqn. (17): 

n= je~' 2+2tt e 2 e 2 dx 
or the alternative summation definition in Eqn. (16) into Eqn. (17): 



(18) 



('-*)' fi 
e 2 e 2 



H a (x) 
ml 



t m dx 



(19) 



or: 
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n =-rSf e 2 e 2 

m! J 



(20) 



Note that the integral in Eqn. (20) is Y m of Eqn. (1). 

Setting the right hand side of Eqn. (18) equal to the right hand side of Eqn. (20) 
produces: 



\e t+2tx e 2 e 2 dx = — Y{e 2 e 2 HJx)dx 
3 ml J 



(21) 



Noting that exponents add, collecting all exponents within the integral on the right hand 
side of Eqn. (17), and setting this to equal E: 



E = - 



t 2 ~ \ ( Z ~ X f X 2 



or: 



E = - 



r-2tx+ 



[z 2 -2zx + x 2 } 



or: 



2 2 2 

E = -(t 2 -2tx + - — ZX + —+-) 
2 2 2 



or: 



E = -{t 2 - x{2t + z) + — + x 2 ] 
2 

Define: 

z = 2y 

and rewrite Eqn. (25) using Eqn. (26): 

E = -[f 2 - 2*(f + y) + 2y 2 + x 2 ] 

Also define: 

F = -[x-(y+t)} 2 

or: 

F = -[x 2 -2x(y + t) + (y+t) 2 ] 

or: 



(22) 



(23) 



(24) 

(25) 

(26) 
(27) 
(28) 
(29) 
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F = -[x 2 -2x(y+t)+y 2 + 2yt + t 2 ~\ (30) 

Thus, 

E-F = -[y 2 -2yt] (31) 
Substituting Eqn. (27) into the left hand side of Eqn. (17) and defining this as R: 

R=\e-^ y ^dx (32) 

or: 

R = e iy2 ~ 2yt) le F dx (33) 

so that: 

R=e^-™\e**^dx (34) 

and, since y and t are constants, thus dx=dy : =0, then R may be expressed as: 

R = e -(y-2^) J e -f-(^)] 2 rf[jc _( y + 1)] (35) 

Setting: 

w = x -(^ + ^ (36) 

then: 

R^e-^^je^dw (37) 

but: 

jV^rfw^ (38) 

Thus, Eqn. (37) becomes: 

R = Je^- 2 ^ (39) 



Substituting Eqn. (39) into Eqn. (20): 

Y^-Y m =R = ^e {y2 - 2y,] 
m 

A Taylor Expansion of e 2yt provides: 



^.. 1+ 3zt + M + _ + fi2i: + _ 

1! 2! w! 



(40) 



(41) 



or: 
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m 




(42) 



Substituting Eqn. (42) into Eqn. (40) yields: 



+m 1 

ml 



m 



ml 



(43) 



Equating coefficients of like terms: 




(44) 



m! 



such that: 



(45) 



Substituting Eqn. (26) in Eqn. (45): 



(46) 



Upon integration of Eqn. (5), a polynomial form, I m (R) m , results and must be 
transformed back into the Hermite polynomial form, J m H m (R) 9 where J m is related to I m 
by a linear transformation which depends on the order, m, of the polynomial. 

Eqn. (5) can be solved exactly. Further, the two sides of Eqn. (5) may be 
represented in the same form. Equating coefficients of like terms enables direct 
relationship of the parameters of the two measurements (actual and captured), which then 
determines the required algorithms. 

The present invention makes use of the observation that most optical and electro- 
optical devices, as well as many other measurement systems, are characterized by a 
Gaussian PSF. Only the constants in the Gaussian vary from system to system. Note that 
the PSF appears within the convolution integral. Three aspects of the solution merge 
here. 

First, Gaussian functions are part of a larger class of functions and are Hermite 
functions of zero order. Second, the actual (undistorted) measurement or image usually 
may be represented by a Fourier-Hermite series. Note also that the actual (undistorted) 
measurement or image is also defined within the convolution integral. Finally, the 
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combination of the first two aspects allows the closed-form, analytic evaluation of the 
nearly general convolution integral without approximation, with the exception that the 
series are finite in length. In the above case, it is an optical convolution integral, but the 
process would work for any device the response of which may be described with a 
Gaussian distribution. Thus, using this approach, the common method of "discretizing" is 
avoided and there is no large system of linear equations that may be incompatible and 
may lead to meaningless results. 

' A preferred embodiment of the present invention, i.e., an evaluation of the 
convolution integral, yields a form proportional to a Gaussian function times a power 
series involving the unknown parameters of the actual (undistorted) measurement or 
image. Using a series of special transformations, the side of the equation representing the 
actual (undistorted) image may be changed into a separate Fourier-Hermite series. 

Since the side of the equation representing the captured (distorted) measurement 
or image may be represented by a Fourier-Hermite series also, the two sides will have 
identical structure. Equating coefficients of like terms now directly relates the parameters 
of the two measurements or images (actual and captured), and permits selection of the 
required algorithms. The only non-numeric constants in the results are the parameters of 
the PSF that have been measured for the individual device or device type. Developing 
this concept further by substituting Eqn. (6) for p(x) into Eqn. (5) yields: 

(z-xf 

D(z)=\e 2 ° I{y)dx (47) 

Substituting Eqns. (3) and (4) into Eqn. (47) yields: 
D(z)=\e >° e 2 £/ m // m (>0^ 

or 

D(z) = P^I m je >° H m (y)e 2 dx (49) 
If the integral on the right side of Eqn. (49) are designated Y m then: 

Y m = \e 2 " e >H m (y)dx (50) 
Since Y m represents a value for the m' H mode, the sum substituted in Eqn. (49) is: 
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To evaluate the integral, choose: 

e 1° =e 2 =e 2 (52) 



Thus, 

2a 2 



2 2 



or 



Let: 



so that 



and let 



then Eqn. (56) may be re-written as: 



(53) 



= - (54) 
a 

Substituting Eqns. (4) and (54) into Eqn. (3): 

l_ i 
/W = ^&A(^ (55) 

and the integral in Eqn. (50) becomes: 

y m = j e 2a e 2aj^ (x)fl -2^ (56 ) 



v = a 2 x (57) 



dx = a 2 dv (58) 



t^a'z (59) 



Y m =a>\e 2 e*H m (v)dv (60) 

Eqn. (60) is now of the form of the convolution integral where a = L The result 
obtained for the convolution integral is: 

I t 

Y m =x*eU m (61) 

and for a ^ 1 is: 
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10 



15 



i Jl 
Y m =(ax)U *t m 

Substituting the definition of r from Eqn. (59) into Eqn. (62) yields: 



Y m =(axye 



1 

4a 



f iY" 

za 2 

) 



To represent as a Fourier-Hermite series, let: 



2 z 2 



2 " 4a 



or 



7? = 



(2a) 



so that Eqn. (63) may be written: 



m R' 



Y m ={nayi 2 e 2 R" 
Multiplying (icaf by P" 1 or (2jrb)'' A yields: 



Y m =(—) 2 l 2 / 2 R m 
m \2b 



Inserting Eqn. (67) in Eqn. (51) produces: 



^)=^] 2 IU2F/'W 



(62) 



(63) 



(64) 



(65) 



(66) 



(67) 



(68) 



This expression now approaches the form of a Fourier-Hermite series because R n 
can be expressed as a linear sum of weighted Hermite polynomials. Rewriting Eqn. (68): 



(69) 



Note that the Fourier-Hermite series can represent an almost arbitrary function in 
the same sense that the common sin(x), cos(x) series of the Fourier Series does. 
20 The Fourier-Hermite representation of D(R) has the general form: 

D{R) = f j D m e~H m (R) (70) 
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where values of D m are constants and the m th value of H m (R) is represented by the m th 
order Hermite polynomial. D m values are determined from H m (R) by: 

D m = (m\2mn)- X J^D(R)e 2 H m (R)dR (71) 

In practice, a finite number of terms in the series of Eqn. (70) is used to represent 
5 D(R). Assuming that finite number to be M, Eqn. (70) may be re-written as: 

D{R) = Y.D m e^H m {R) (72) 

The number of terms, M, selected here will also determine the number of terms 
for the representation of the actual (undistorted) measurement or image, I(x). The value 
chosen for M determines the specific structure of the algorithm that the procedure 
10 produces. This provides some flexibility for the user in meeting the possibly competing 
requirements of computational efficiency and image "sharpness," for example. 

Based on the "basic" convolution equation of Eqn. (5), substitute the RHS of Eqn. 
(69) for the LHS of Eqn. (5) and the RHS of Eqn. (72) for the RHS of Eqn. (5), yielding: 

M ( n ^2 Jl- M 2 

e 2 ZAA(*H ^ 6 2 I 2 J » H »W (73) 

m=o \£DJ m=0 

15 Because both sides of Eqn. (73) now have identical structures, coefficients of like 

terms may be equated. A set of linear equations relating parameters of the captured 
(distorted) and actual (undistorted) images may be produced. Using the set of Eqns. (11) 
and (12) provides a system of linear equations that can be solved for J m and then 7 m , 
yielding the exact solution to the deconvolution. Using the common multiplier, E, and 

20 Eqn. (14), Eqn. (15) may be expressed alternatively as: 

4=l4A=£X4A (74) 

m=0 m=0 

representing the solution of the deconvolution, where A m is the "structural" part of the 
relationship between I m and J m . Note that I m represents the coefficients of the actual 
(undistorted) measurement or image and J m represents a set of dummy variables. It is 
25 desirable to be able to move back and forth between the two sets. 

Refer to Fig. 2, a flow chart of a process yielding accurate representations of 
actual measurement or image data by deconvolving data collected by systems, such as 
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optical detectors, whose point spread function response follows a Gaussian distribution. 

The process comprises: 

establishing 201 a mathematical relationship as a general convolution integral 
equating the actual measurement or image data to the collected measurement or 
image data; 

selecting 202 a Fourier-Hermite function; 

employing 203 a generating function for Hermite polynomials for expanding a 
Fourier-Hermite series of the Fourier-Hermite function to establish a linear 
mathematical relationship between the actual and the collected measurement or 
image data; 

selecting 204 an order, m, of the Fourier-Hermite polynomial such that m also 
determines the number of terms for the representation of the measurement or 
images; 

expanding 205 the actual measurement or image data in a Fourier-Hermite form 

with unknown (to-be-determined) coefficients by employing a series of special 

transformations to convert the side of the mathematical relationship representing 

the actual measurement or image to a Fourier-Hermite series; 

expanding 206 the collected measurement or image data in a Fourier-Hermite 

form; 

equating 207 coefficients of like terms on each side of the mathematical 
relationship to relate the coefficients of the actual and collected measurement or 
images; 

selecting 208 a linear convolution algorithm from which a general convolutional 
equation is derived; 

solving 209 the resultant general convolution equation exactly, in closed 
analytical form, as a linear convolution equation, such that using this process to 
evaluate the general convolution integral yields a form proportional to a Gaussian 
function times a power series that is defined with a finite number of terms, m, and 
incorporates the unknown coefficients of the actual measurement or image data as 
presented in a Hermite function, such that a solution in closed analytic form 
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provides a satisfactory solution to the general convolution integral without 
approximation; and 

performing 210 an analytic decorivolution of the resultant linear convolution 

algorithm by employing classical deconvolution algorithms solely; 

evaluating 211 the adequacy of the result such that the analytic deconvolution 

yields a solution having acceptable error levels; 

if adequate, stop 212; 

if desirable to achieve greater accuracy, input 213 to an iterative deconvolution 
process. 

As with other embodiments, this process may be used to initiate any conventional 
iterative nonlinear deconvolution technique to further refine the accurate representations. 
Note that fewer iterations are needed because of the accurate starting point provided. 



For simplicity, a process is described for one-dimensional (1-D) imagery 
collection such as produced by the HYDICE sensor. With these types of sensors, a single 
line of optical detectors/pixels oriented transverse to the direction of platform travel scans 
the angular line of sight electronically in a "whiskbroom" operation. McKeown, D. M. et 
al, Fusion of HYDICE Hyperspectral Data with Panchromatic Imagery for Cartographic 
Feature Extraction, IEEE Trans. Geoscience and Remote Sensing, 37, No. 3, pp. 1261- 
1277, May 1999. The individual 1-D rows are later combined in the HYDICE imaging 
process to form a matrix for building 2-D images. 



EXAMPLE 



SAMPLE CALCULATION 



Assume Jo, Ju and J 2 represent a set of dummy variables and 



A(R) = 2X# m (JO = J 0 H 0 + J X H X + J 2 H 2 



(75) 



The Hermite polynomials of orders one, two and three are: 



//„(*) = 1 
H t (R) = 2R 
H 2 {R) = 4R 2 -2 



(76) 



By substitution: 



B(R) = = I i m R m = (Jo -2J 2 )+(2J l )R+WR 2 



(77) 
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It follows that values for/ w are given by: 

I 0 = (J Q -2J 2 ) 



and conversely, 



A -2/, 
/ 2 =4J 2 



1/ 

•'o — A) + 2 2 
' 2 1 



(78) 



(79) 



Dealing with Eqn. (73) again: 



m=0 



m=0 



(80) 



or 



^ „ \ 



I i D a H m (R) = C^J m (2yH n (R)= ±\ %J m (2y H m (R) 



m=0 m=0 

In expanded form, Eqn. (81) reads: 
D 0 H 0 (R) + D l H l (R) + D 2 H 2 (R) = C 



J 0 H 0 (R) + J X (2YH S {R) + 2J 2 H 2 (R) 



Equating coefficients of like terms involving Hermite polynomials: 



D 0 = C/ 0 



(81) 



(82) 



A=C/,(2) 2 
D 2 =2CJ 2 



(83) 



To obtain the solution of the convolution problem, substitute Eqns. (78) into Eqns. (82): 



D 0 = C 



I n +-l 



(84) 



A = 
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To solve the deconvolution problem, these relations are inverted: 





|A 

_c _ 




A 
. c. 




/,=(2)^ 


"Al 
c 




(85) 


h = 


2 


D 
C 


2 
i 









where, remembering Eqn.(6): 

a is a measure of the relative width, and 

b is an inverse measure of the relative height of the Gaussian system PSF and 
Do, Dj, and A? are the full set of three actual (undistorted) image parameter 
measures- 
While the invention has been described in terms of its preferred embodiments, 
those skilled in the art will recognize that the invention can be practiced with 
modifications within the spirit and scope of the appended claims. For example, although 
the system is described in specific examples for topography, it will operate in areas of 
medicine, communications, and even business models where deconvolution is a preferred 
procedure. Thus, it is intended that all matter contained in the foregoing description or 
shown in the accompanying drawings shall be interpreted as illustrative rather than 
limiting, and the invention should be defined only in accordance with the following 
claims and their equivalents. 
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